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Abstract. 

Theories of gravity other than general relativity (GR) can explain the observed 
cosmic acceleration without a cosmological constant. One such class of theories of 
gravity is f(R). Metric f(R) theories have been proven to be equivalent to Brans-Dickc 
(BD) scalar-tensor gravity without a kinetic term (oj = 0). Using this equivalence and 
a 3+1 decomposition of the theory it has been shown that metric f(R) gravity admits 
a well-posed initial value problem. However, it has not been proven that the 3+1 
evolution equations of metric f(R) gravity preserve the (hamiltonian and momentum) 
constraints. In this paper we show that this is indeed the case. In addition, we show 
that the mathematical form of the constraint propagation equations in BD-equilavent 
f(R) gravity and in f{R) gravity in both the Jordan and Einstein frames, is exactly 
the same as in the standard ADM 3+1 decomposition of GR. Finally, we point out 
that current numerical relativity codes can incorporate the 3+1 evolution equations of 
metric f(R) gravity by modifying the stress-energy tensor and adding an additional 
scalar field evolution equation. We hope that this work will serve as a starting point 
for relativists to develop fully dynamical codes for valid f(R) models. 



PACS numbers: 04.25.D-,04.25.dk,04.30.-w,04.50.Kd 
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1. Introduction 

Since its formulation, Einstein's general theory of relativity (GR) has withstood 
extensive experimental and observational scrutiny using tests that range from millimeter 
to solar system scales (see [1] and references therein). The discovery of the late-time 
cosmic acceleration j2j [3] was a surprise, but one which could be modelled within the 
minimally extended framework of ACDM [H [5] - GR with a positive cosmological 
constant. To this day this simple model remains in very good agreement with data 
from all competitive probes El E], which imply that approximately 70% of the 
energy density of the universe is made up of a component which does not cluster and 
has an equation of state with pressure approximately equal to minus the energy density. 
While the simplest model for this component is indeed the cosmological constant, from 
the point of view of particle physics, its value implied by the measurements of the 
cosmological expansion is extremely low and requires a very high level of fine tuning. 

A number of alternative models for dark energy have been proposed, most of which 
suffer from a similar fine-tuning problem to ACDM (see the review [10]), but at least 
provide a set of alternatives against which to test the ACDM hypothesis. In this spirit, 
it is possible to imagine that, rather than proposing the existence of a new, exotic form 
of energy density, it is the theory of gravity which we use to interpret the cosmological 
data that must be modified. 

There are a number of proposed gravity theories which modify the dynamics at large 
distances, and metric f{R) theories of gravity (see, e.g., [HI [12] and references therein) 
comprise one such class of modifications to GR. This class has attracted considerable 
attention in recent years, perhaps due to the simplicity of the modifications. Further 
motivation for the study of f{R) gravity is reviewed in [12J; for other interesting 
alternatives, see [11] for Gauss-Bonet gravity, [13] for conformal gravity and [H] for 
Brane- World gravity. 

The f(R) formulation arises from a simple replacement of the Ricci scalar (R) in 
the Einstein-Hilbert action, 

S= ^j ^ d ^ R ~ 2A ) + S m{9^m): (1) 

where g is the determinant of the metric tensor g^, A the cosmological constant, S m the 
matter term in the action, and ip m collectively denotes the matter fields, by an arbitrary 
function of the Ricci scalar, i.e., 

S = 7^- / y/^gd^xfiR) + Smig^ipm). (2) 
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Note that throughout this work we adopt geometrized units, where G — c — 1. 

From Equations ([1]) and (j2J) it is clear that GR is recovered for f(R) = R — 2 A. 
In metric f(R) theories the connection symbols are chosen to be the Christoffel 

symbols associated with the metric tensor, so that the action is a function of only the 
metric tensor and its derivatives. As a result, in metric f(R) gravity only the metric 
tensor is truly dynamical. In Palatini f(R) gravity the connections are considered 
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independent of the metric tensor, so that the action is a function of both the metric 
tensor and the connection symbols. Thus, in Palatini f(R) both g^ v and ^r*^ are 
dynamical fields (see also [15] for a new class of models which interpolate between the 
metric and Palatini formulations). In this work we are concerned with metric f(R) 
gravity only. 

Early work on f(R) theories [161 El HI EEH] was mainly concerned with high- 
energy corrections to general relativity and their influence on the early universe (see in 
particular [19] where the first f(R) model of inflation was proposed). The discovery 
of cosmic acceleration [21 [3] renewed the interest in f(R) models, but now with 
modifications in the infra-red. A number of alternative models to GR have been 
proposed [201 I2H 1221 [231 [21] ■ However, it was later shown that these models neither 
satisfy local gravity constraints [25l [26j [27] nor give rise to a standard matter-dominated 
era [2E1E9]. 

General conditions for the cosmological viability of f(R) models were derived in 
[30] and it was later realized that the so-called Chameleon mechanism - the scalar degree 
of freedom becomes massive in dense environments and light in diffuse ones - can allow 
f(R) gravity to satisfy Solar-System constraints [2H [32] . The key consequence of the 
Chameleon mechanism is that the modification to the metric inside galactic haloes is 
suppressed: gravity returns to its general-relativistic behaviour. The functioning of 
the Chameleon mechanism has also been confirmed via N-body simulations of large- 
scale cosmological structure formation in [331 Ell EHl [361 EI] > where it was shown that 
predictions for cluster abundance and the matter power spectrum return at small scales 
to those calculated within the ACDM framework. 

A number of models that satisfy both Solar-System and cosmological constraints 
have been proposed in [311 [321 EH EH HOI HH H21 S3], and it is now known that for an 
f(R) theory to be viable the following four constraints must be met [TT] : 

(i) f,R > for R > R , where i? is the cosmological value of the Ricci scalar today. 
This condition is necessary for guaranteeing that the new scalar degree of freedom 
is not a ghost - a field with negative kinetic energy. 

(ii) f(R) — > R for R ^> Rq. This condition is necessary for the presence of a matter- 
dominated era and to evade solar-system constraints. 

(iii) f,RR > for R > Rq in the presence of external matter. This condition ensures 
that the matter-dominated era is the stable solution for cosmology and that the 
solutions which satisfy solar system constraints are stable. 

(iv) < Rf,RR I f,R.\ r =-2, where r = —Rf,R/f. This condition is necessary for the 
stability and presence of a late time de Sitter solution. 

The existence of these requirements is a result of the fact that in f(R) gravity 
the Ricci scalar is a full dynamical degree of freedom, which must behave in a manner 
similar to the Ricci scalar in GR, where it is controlled through a constraint (R = —87rT). 
These conditions ensure that in high-density environments, the so-called high-curvature 
solutions, where R ~ Srrp, are stable. 
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An additional constraint that any theory of gravity must satisfy is the existence 
of stable relativistic (neutron) stars. It was originally pointed out in |44j . that many 
models of f(R) theories reach a curvature singularity at a finite value of the scalar degree 
of freedom which is not protected by the existence of a potential barrier. This value 
of the scalar field may be attained in the presence of relativistic matter. This same 
idea was used in [15] to argue that it is not possible to build spherically symmetric, 
i.e., non- rotating, relativistic stars in f(R) theories of gravity. These works stimulated 
further interest and eventually numerical models of spherical relativistic stars in f(R) 
gravity were explicitly constructed in [l6l H7J HH]. There it was shown that building 
numerical models of neutron stars in f(R) gravity is very sensitive to the treatment of 
boundary conditions. 

To our knowledge a stability analysis of non-rotating equilibrium models of neutron 
stars in the context of f(R) theories has not been carried out yet. One may expect 
that the stability properties of relativistic stars in viable f(R) gravity are the same 
as those in GR, because of condition 3 above. However, given the subtleties that 
arise in obtaining relativistic stellar configurations in f(R) theories due to the effective 
scalar degree of freedom it is natural to expect that the back-reaction of the scalar 
field will affect the stability, too. In addition, it would be interesting to explore the 
existence and stability of rotating neutron stars and how f(R) gravity affects the 
criterion for the onset of the bar mode, r-mode and other non-axisymmetric instabilities 
[4~9| \50\ I5TI |52| [53| 154"] . Furthermore, it is intriguing to study gravitational radiation 
arising from compact stars, both in isolation and in binary systems. Included in this 
list are neutron star - neutron star [55], black hole-black hole [56], black hole-neutron 
star ^3 EEl IBOl ISH 15H1 15H IEEI EZl IES1 EE EH IZ21 IZ3 ISl IS] and. white- 
dwarf-neutron star binaries [761 [77J . 

Some of these studies can be carried out analytically via perturbation theory, 
and some require direct numerical simulations. One of the main points we make in 
this work is that current numerical relativity techniques (see texts by Baumgarte and 
Shapiro [78] and Alcubierre [79] and references therein), i.e., the solution of the Einstein 
equations by computational means, should be able to handle the equations of f(R) 
gravity straightforwardly. In particular, the minimum requirement is to modify the 
stress-energy tensor and add a new scalar field evolution equation. However, to achieve 
long-term stable numerical integration of any set of partial differential equations, well- 
posedness of the Cauchy (or initial value) problem must be guaranteed. 

Unlike GR, the field equations of metric f(R) gravity in the so-called Jordan frame 
are 4th order (see Section [2]). Nevertheless these theories can be cast in 2nd-order form, 
by promoting (the derivative of f(R) with respect to R) to an effective dynamical 
scalar degree of freedom. Alternatively, metric f(R) gravity can be reduced to second- 
order form by a transformation of the f(R) action to a Brans-Dicke (BD) [51)] action 
with lo = [25J. This means that f(R) gravity is equivalent to BD gravity without a 
kinetic term. Exploiting this equivalence and the 3+1 decomposition approach of [81J, 
it was demonstrated in [82] that metric f(R) gravity admits a well-posed initial value 
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problem. As in 3+1 GR, to solve the initial value problem, first one solves the 3+1 
constraint equations to obtain initial data and then uses the 3+1 evolution equations 
to advance the initial data in time. For this approach to yield a consistent solution 
of the covariant (4D) field equations, the 3+1 evolution equations must preserve the 
constraints of the theory. To prove this one has to derive the evolution equations of the 
constraints, which are often referred to as the constraint propagation equations, and 
show that if the constraint equations are initially satisfied, they must be satisfied for 
all times. To our knowledge this has never been demonstrated for a 3+1 formulation of 
f{R) gravity and in this work we show that this is indeed the case. 

To date there are two methods for deriving 3+1 constraint propagation equations. 
One approach is to take the time derivative of the constraint equations in 3+1 form 
and then replace the time derivatives of all dynamical variables by using the evolution 
equations for these variables. We call this the 3+1 or "brute force" method. This is 
a rather tedious approach and to our knowledge, it has been performed in GR only 
for vacuum spacetimes in [83J. A pedagogical example that explains the "brute force" 
method is given in section II of [83], and more involved applications involving Maxwell's 
equations can be found in [S3]. 

The other approach, which is more elegant, takes advantage of the Bianchi 
identities. We call this the Frittelli method [85] (see also |79j). 

However, the equations derived in [85] were not cast in pure 3+1 language. Here 
and throughout this paper by "pure 3+1 language" we mean that a given equation 
is written solely in terms of scalars and purely spatial objects and their derivatives. 
Yoneda and Shinkai [861 [87] have derived the Arnowitt, Deser, Misner (ADM) constraint 
propagation equations in pure 3+1 language but they did not indicate how they arrived 
at their result. 

In this work we employ the Frittelli approach to derive the constraint propagation 
equations of f(R) gravity and cast the resulting equations in pure 3+1 language. We 
show that the mathematical form of the constraint propagation equations is the same 
as that of the standard ADM formulation of GR. We also demonstrate that this result 
holds true both in the Jordan and the Einstein frames of metric f(R) gravity, as well as 
for the BD-equivalent version of metric f{R) gravity. Finally, we compare our equations 
with published results of the constraint propagation equations derived using the 3+1 
approach and show that the expressions obtained via both approaches agree. 

While none of our results are surprising they serve to prove that f(R) gravity is self- 
consistent. Moreover, it is revealing to demonstrate how previous results from GR can 
be extended to alternative theories of gravity and the consistency between alternative 
approaches. Finally, obtaining the extended constraint propagation equations in pure 
3+1 form may prove useful for performing 3+1 numerical simulations, where constraint 
preservation can be used as a check on the integration. 

This paper is organized as follows. In Section [2] we review the field equations of 
generic metric f{R) models. In Section [3] we provide a simple pedagogical argument 
(see also [88j[78]) to demonstrate the basic idea of constraint preservation in the context 
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of GR. In Section H] we review the 3+1 decomposition of the BD-equivalent metric f(R) 
equations. In Section[5]we employ the Frittelli method and use the results of Section0]to 
derive the 3+1 metric f(R) constraint propagation equations. In Section [6] we cast our 
generalized evolution equations of the constraints in pure 3+1 language. In Section [7] we 
argue that the 3+1 constraint propagation equations of f(R) gravity in both the Jordan 
and the Einstein frames can be cast in the same form as that in the 3+1 BD-equivalent 
version of f(R) theories. Finally, we summarize our work in Section [HJ 



2. f(R) field equations 

As in GR, the fundamental quantity in f(R) gravity is the spacetime metric tensor g a p 

ds 2 = g a pdx a dx 13 , (3) 

where ds is the line element, and x a denote the spacetime coordinates. Here and 
throughout this paper Greek indices run from to 3, while Latin indices run from 
1 to 3. 

The goal of the theory is to determine the metric given a mass-energy distribution. 
Because of the existence of an additional scalar degree of freedom in the gravitational 
field sector, it is possible to formulate the field equations of f(R) theory in many ways, 
depending on the amount of mixing between these two fields. We will discuss three such 
formulations: the Jordan frame, the Einstein frame, and the BD-equivalent formulation. 

The Jordan frame and Einstein frame formulations have different metrics as 
dynamical variables. The two metric tensors are related via a conformal transformation 

g^ = Vt 2 g^ u (4) 

where Q is the conformal factor, and g^ here denotes the metric in the Jordan frame. 
Note that Equation (jlj) is equivalent to a transformation of units |89j . 

In this section we review the field equations of f(R) gravity in both the Jordan and 
Einstein frames, as well as those of the BD-equivalent form of f(R) gravity^. 



2.1. Jordan Frame 

The action (j2J) is called the Jordan frame action. An action is said to be in the Jordan 
frame, if the dynamical metric tensor in the action is the metric whose geodesies particles 
follow, i.e, the physical metric. The Jordan frame is the one in which the definition of 
the matter stress-energy tensor is 

where 8S m /b~g^ u is the functional derivative of S m with respect to g^ u . For example, it 
is in this frame that the stress-energy tensor of a perfect fluid has the form 

= (P + P)u,u v + P 9lun (6) 
| For various Hamiltonian formulations of f(R) gravity see |90| . 
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where p is the total energy density of the fluid, P the fluid pressure, and the fluid 
four velocity. 

Varying the action (j2J) with respect to the metric yields the f(R) field equations 
[TT| [12] in the Jordan frame 

^ = 8*7*?, (7) 

where T^r is the matter stress-energy tensor and 

5V = Fi^ - ~/<? M „ - V M V,F + g^DF, (8) 

and where F = f,R. Note that for brevity we have dropped the argument of both f(R) 
and F(R). Clearly GR is recovered for f = R — 2 A, in which case Equations (J7|) and 
(EJ) yield 

£^ = G Mi/ + A^ = 87rT^), (9) 

where G MV is the Einstein tensor. 

Equation (0) is 4th-order due to the term V^V^F. However, if we take the trace 
of Equation (JTJ), we obtain 

3DF + FR-2f = 8vrT (m) , (10) 

where T( m ) = and 

□f = -L^v^rw (ii) 

Equation (ITUj) can be used to promote F(R) into an effective dynamical scalar degree of 
freedom (often referred to as "scalaron"), thus recasting the theory in 2nd-order form. 
Equation (JTJ) can also be written in the following form [39] 

G^ = 8ir(T^+T$), (12) 
where Tpi, can be thought of as a "dark energy" stress-energy tensor, given by 

87rT#>= l - g ^f-R) + v,V v F 

-g^nF+{l-F)R tlv , (13) 

This form of the field equations of the theory is interesting because the Bianchi 
identities VG^ = together with VT^ = 0, imply that 

VTW = 0, (14) 
i.e., the dark energy tensor Tpi is conserved. 
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2.2. Einstein frame 

To obtain the Einstein frame action of f(R) gravity, i.e., an action linear in a Ricci scalar 
R associated with a metric g^, all we have to do is perform a conformal transformation 
on the metric 

= Fg^, (15) 

i.e., the conformal factor Q in Equation (j3J) is Q 2 = F. For the transformation to be 
physical F must satisfy F > 0. Note that this condition is in accord with the first 
condition for cosmological viability of f(R) gravity listed in Section [TJ 
If we introduce a new field 6 such that 



3 InF, (16) 
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then the f(R) Jordan action transforms to [TT 



where 



Se= 1^kG J dAx ^~9 R + S 4> + S ™( F 1 ( ( f ) )9^,^m), (17) 

S^ = J d 4 x^[ - l -g^d^d v <P - V(<f>)\ (18) 
is the scalar field term in the action, and where the scalar field potential is defined as 

= < 19) 

The dynamical metric tensor in the Einstein frame is not the physical (g^ u ) but 
the conformal one {g^). However, the matter still follows the geodesies of the physical 
(Jordan) metric. Variation of the matter action with respect to g^ u yields 

f(m) = 2 = J_ T (m) / 9f) x 

which is no longer independent of the scalar field <fi. 

Variation of the action (fTTj) with respect to yields the scalar field equation 

<"> = (), (21) 

where T (m ) = g^fff and 

□0 = -^d^^l~g^d v <P). (22) 

Equation (121]) implies that the scalar field is directly coupled to matter. 
Finally, variation of the action (jTTj) with respect to g^ yields 

^ = 87r(f(™)+fW), (23) 

where the scalar field stress-energy tensor is 

1 



W = d„4>d u (t> - ~ QiiV 



2 g a$ d a (/>dp(f> + V(<f>)\. (24) 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 9 
Note that in the Einstein frame V^TJu™ 7^ 0; instead we have 

v^ = v*(fM + rW) = o, (25) 

where V M is the covariant derivative associated with g^ v . It can also be shown that [11] 
V"f#> = ~^fV u <t>, V"f = ^fV v <t>. (26) 

2.3. Equivalence with Brans-Dicke gravity 

Another way to cast f{R) gravity into second-order form is to express the theory as a 
BD theory. To show that f(R) gravity is equivalent to BD gravity with a potential, the 
following action was considered in [25J: 



S = — / V^d'xifix) + f, x ( X )(R ~ X)\ + S m . (27) 

107T J 

Varying the action with respect to x yields 

f, xx ( X )(R-x) = 0- (28) 

Thus, if f, xx (x) 7^ (in agreement with condition 3 in Section [TJ, then 

X = R. (29) 

Hence, Equation (1271) recovers the Jordan frame f(R) action (j2J). If we now let 
= fix (x)> Equation ( 1271) can be written as follows 

S = JL / ^ 4 X U - 1/(0)1 + S m , (30) 
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where the potential is given by 

V(<f>)= X ((f>)(f>-f(x(<f>))- (31) 

Action (1301) is the same as the original BD action with a potential and without the 
kinetic term (u/2)g' 1L/ d^<j)d u <j), i.e., the BD parameter is u = 0. Varying the action (1301) 
with respect to the metric yields the BD-equivalent f(R) field equations [12] 

8tt 

where 



G^ = -^(T^+rW), (32) 



SttTW = V„V„0 - ^(D0 + ^(0)) (33) 
Taking the trace of Equation (|32|) we can replace i? in Equation ( 129]) to obtain 

3D0 + 2VU) - (j)^- = 8vrT. (34) 
dcp 

Equations (|32|) and (134")) are the BD-equivalent f(R) field equations. 
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3. Constraint Preservation in a spacetime context 

In this section we review the concept of constraint preservation using the standard 
Einstein equations in 4D covariant form, i.e., we do not invoke machinery of the 3+1 
decomposition of spacetime. The reason for doing so is that so far we have written the 
most popular representations of metric f(R) field equations in a GR-like form 

= 8nT^, (35) 

where is an "effective" stress-energy tensor that is conserved, i.e., V M T Ml/ = 0. Thus, 
it is instructive to first consider the Einstein equations in their familiar 4D covariant 
form. 

For the Einstein equations with a cosmological constant we have T^ v = — 
Ag^ u /8n. Since the Einstein equations are second-order partial differential equations, 
the evolution of the 4- metric g a p in time can be determined by specifying g a p and d t g a p, 
everywhere on a three-dimensional spacelike hypersurface that corresponds to a given 
initial time t. Equation ( l35l) can provide us with expressions for d^g a p, which we can 
use to evolve the metric in time. There are 10 metric components and there are 10 field 
equations in (!35| . Hence, it appears that we have the exact number of equations for the 
10 degrees of freedom of the metric. However, the Bianchi identities V / gG a/3 = 0, give 

d t G a0 = -d t G ai - G^ (4) r% - G Q ^ 4 )r%, (36) 

where we set d t = do and where ^T ^ are the Christoffel symbols associated with g a p. 
Since no term on the right-hand-side of Equation ( 136]) contains third time derivatives or 
higher, the four quantities G a0 cannot contain second time derivatives. Thus, the four 
equations 

G^ = 87r7> (37) 

do not provide any information on the dynamical evolution of the metric. They are 
instead a set of constraints that g a p and dtg a p have to satisfy. The only truly dynamical 
equations are the six remaining equations 

Gn = 8tt% (38) 

The apparent mismatch between the number of metric components and the number of 
evolution equations is immediately resolved once we invoke the coordinate freedom of 
GR. The theory is four-dimensional, and hence we can always choose four conditions 
to specify a coordinate system. For example, we can choose the four g p components 
and assign them certain values, or demand that they satisfy a given set of four partial 
differential equations. This way we are left with six independent metric components, 
for which we have the exact number of evolution equations ( 138]) . 

However, solving Equation ( 138]) does not guarantee that the full set of the Einstein 
equations (153]) will be satisfied. For that to be true, Equation ( 137]) has to be satisfied 
for all times. In other words, if one solves Equation ( 138]) starting with initial data that 
satisfy Equation ( l3~Tj) . one has to prove that the constraints are preserved. 
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To demonstrate that this is indeed the case we make use of the Bianchi identities 
in the following form: 

= 0, (39) 

or 

d t s a0 = -di£ ai - £^ (4) rv - £ a ^ (4) rv, (40) 

where 

£ aP = Qa/3 _ 87r f af) . (41) 

If we let C = £ 00 and C % = £ t0 , Equation ( 1401) can be rewritten as 
d t C= -d i C l -C{2^T\ + ^T i 0i ) 

-^(3Wr% + ( 4 )r%-) (42) 
d t c j = - c (4) r j 00 - 2C i(4) r^ - c i(4) rV (43) 

where we have used Equation ( l38l) . £ l i = 0, to obtain the result. Thus, if the constraints 
are initially satisfied, then C = C % = initially and from Equation ( 1421) the time 
derivative of the constraints will be zero and hence the constraints will remain zero for 
all times. Since this conclusion resulted from setting = 0, the previous statement is 
equivalent to saying that the evolution equations preserve the constraints, a result that 
is well-known. 



4. 3+1 Decomposition of f(R) gravity 

Well-posedness of the Cauchy problem in metric f(R) gravity has been demonstrated 
in [82] using the BD-equivalent f(R) formulation. In this section we focus on the BD 
version of f(R) gravity and review the salient features of its 3+1 decomposition that 
will be useful in our proof of constraint preservation. 

The form of the field equations of the theory is that of Equation (1551) . where 

r^ = l(TM + T^). (44) 

The 3+1 decomposition of spacetime is a decomposition of spacetime into space 
and time. To do this, one assumes that the four- dimensional spacetime manifold can 
be foliated by a one-parameter family of nonintersecting spacelike hypersurfaces. The 
parameter of this family of hypersurfaces is taken to be the coordinate time. The 
spacetime metric is then rewritten as [9TJ 

ds 2 = -a 2 dt 2 + j ij (dx i + ftdt){dx j + f3 j dt), (45) 

where a is the lapse function, (3 l is the shift vector, and is the 3-metric on the 
spacelike hypersurfaces, induced by g a p. The lapse function and the shift vector are 
gauge quantities; they dictate how to build the coordinate system and can be freely 
specified. The relation between 7^ and g a p is 

1* p = + rfnp, (46) 
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where 7°^ = g^Jnp, b~ a p is the Kronecker delta, and n a is the future directed timelike 
unit vector normal to the t = const, hypersurfaces. The tensor 7°^ is the operator that 
projects tensors onto spacelike hypersurfaces. 

The field equations can then be decomposed into a set of evolution equations and 
a set of constraint equations by using 7 a /3 and n a . 

Projecting Equation (135]) twice with the projection operator yields the evolution 
equations 

= (G a p - 87rf^) 7 %7^ = G a ^«lt ~ 8ttS^ = 0, (47) 



where 



and where 



ee T Q/37 % 7 ^ = + Sj$ , (4? 



^ - Z^tVA, « } - 7$7VyV (49) 
Using Equation ( 133]) . we can write in Equation ( j4"9j) as follows 

3$ = ^i D ^^ - + ^(0)1, (50) 

where _D M is the covariant derivative associated with 7^. 

Contracting Equation fl35|) twice with n a yields the Hamiltonian constraint 

H ee (G Q/9 - 8nf a/3 )n a n^ = G a pn a n? - 8vrp = 0, (51) 

where 

p^T^rfn^p + pW, (52) 

and where 

p ee rJ'nV, p W ee T^V (53) 
Using Equation (1331) we also obtain 

p W = -L^n'V M V^ + (D0 + V(0))]. (54) 

Contracting Equation (1351) once with n a and projecting once with j a p yields the 
momentum constraints 



M„= - (G afi - 87rT a ^)n a 7^ 

= - G^n Q 7 ^ - 8tt^ - torjjf> = 0, ( 55 ) 

where 



J, = -iJW,, J'f = -3gW„ (56) 
and where from Equation (I3"3"j) we find 

# } = ^>VmV«V^. (57) 
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We can now write Equation (1351) as a linear combination of the evolution and the 
constraint equations. 

G«f> _ Snf a ? = (G^ - 8nf^)5 a ^ u 

= (G^- SttTH ( 7 % - n a n,) tf„ - n?n u ) (58) 

where we used Equation ( HST) in the second line to replace the Kronecker deltas. By use 
of Equations (T4"Tj) . fl5T|) and fl55|) . Equation (158]) becomes 

- 8vrf a/3 = E aP + 2n (a M /3) + ffn a /. (59) 

This last equation has been derived by Frittelli, (see Eq. (9) in [85]) for GR. Here we 
have shown that this equation is valid in f(R) gravity, too, provided that appropriate 
definitions of H and M % are given. 

As was shown in [85] , setting E al3 = yields the evolution equations of the original 
ADM formulation [91], whereas setting E al3 = 7 a/3 if yields the evolution equations of 
the standard ADM formulation J92J [T9J, [78]. Following the parametrization of [85] we 
set E a/3 = A7 a/3 if , so that in the ADM language A = corresponds to the original ADM 
formulation, while A = 1 to the standard ADM formulation, except that here we deal 
with f(R) 3+1 formulations. 

It is now evident from Equation ( 1591) that if M a = H = and E a/3 = A7 a/3 if , then 
the f(R) equations are satisfied. 

By introducing the extrinsic curvature K~ij 

Kij = —^£njij (60) 

where £ n stands for the Lie derivative along the timelike unit vector n a , using the 
Gauss, Godazzi and Ricci equations (see e.g. Equations (2.68), (2.73), (2.82) in [78]). 
and adopting the usual coordinate basis where 

= {a' 1 , -a- l f3 l ) (61) 

one can derive the evolution and constraint equations in 3+1 form, which (for A = 0) 
are presented in [82] and we do not repeat them here. 

A subtlety that must be addressed for our purpose and which is pointed out in 
[8~T| [82] is that to remove the time derivatives of the scalar field from the sources 

pMjjP one introduces the gradients of <fi as new dynamical variables 

n = £ n <p = n^0, (62) 

Qn = D^. (63) 

The □</> operator in the sources SjiJ , pMjfi can be removed by use of Equation 
( )34|) . Furthermore, Equation (13^|) in combination with Equation (|62|) . which can be 
written as 

aU = d t <f) - FQt, (64) 
can be used to derive the evolution equation for II. Eventually, one finds [8T] 

£ n n = UK + CyDiQn a) + A<7 - U(j>. (65) 
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where K = j^Ky. 

To promote Qi to a dynamical variable we take a time derivative of Qi and using 
Equation (IMj) we obtain 

d t Qi = £pQi + Di(an), (66) 

where 

£pQi = P s d s Q l + Q s dif3 s . (67) 

The introduction of new variables Qi, introduces an extra constraint, which the 
evolution equations have to satisfy 

Ci = Qi- D i( p = 0. (68) 

In addition to this, the ordering constraint 

Cy = DiQj - DjQ, = 0. (69) 

has to be satisfied, too. 

Thus, constraint preservation means that the evolution equations must preserve all 
the constraints of the 3+1 decomposition, i.e., Equations ( |5T|) . ( 1551) . ( 1681) . and ( |69f) . 



5. Constraint Propagation Equations of 3+1 f(R) gravity 

The backbone of the Frittelli approach is to express the field equations in the form 
of Equation ( |59i) and plug it in the Bianchi identities in order to derive the evolution 
equations for the constraints, assuming the evolution equations are satisfied E^ u = 
X^^H . So far we have extended the Frittelli approach to general metric f{R) gravity. 
Since the form of Equations ( |59l) is the same as in [85], the derivation of the 3+1 BD- 
equivalent f{R) constraint propagation equations is precisely the same as that in [85J, 
which is valid for GR, and to which we refer the interested reader for more details. Here 
we only sketch the derivation and write the result. 
The Bianchi identities are 

V M (G^ - 8vrf = 0. (70) 

or, equivalently, after substituting Equation (159]) in Equation (jTOjl 

V M (£^ + 2n { ^M u) + Hn^rf) = 0. (71) 

To find the evolution of the Hamiltonian constraint we contract Equation (17T1) with 
n a and after some algebra we find 



= - E^DyUfj, - 2n u M^V u n fl - D V M' 
- n u V u H — HD u n u . 



(72) 
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To find the evolution of the Momentum constraint we project Equation (I7ip with 
7°« and after some algebra we find T 



= D^E*"* + nTE^V^ + ^V M M Q - n a M u n^V^n v 

+ M a D fl n fl + M^D^n® + Hn^V^n 01 . (73) 

Proof that our equations are correct will be provided below when we cast the 
constraint propagation equations in pure 3+1 language and compare our result with 
published results in the literature obtained via the 3+1 approach. 

Using the following identities 

-fHDpn,, = HDX, (74) 
n^HV^nu = Hn^V^n 01 , (75) 

and substituting E^ u = Xj^H in Equations (172]) and (!73|) we find that the evolution of 
the constraints is given by 

n^V^H = -2n»M u V >„ - D^M" - (1 + X)HD fl n^ } (76) 

n^V^ikT = - + n u M a n p V p n a - M v D i jv a 

- M^D^n u - (1 + X)Hn^V fl n u . (77) 

These last two equations have the same mathematical form (except for a factor 
of 2; see footnote in page [15]) as those derived in [85] that applied to the case of GR, 
i.e. f(R) = R. Here we have proven that the form of the hamiltonian and momentum 
constraint propagation equations is the same for both vacuum (T^ = 0) and non- 
vacuum spacetimes (T^ ^ 0), and that it is independent of the f(R) function, because 
we have absorbed all terms that depend on these quantities in the definition of the 
hamiltonian and momentum constraints (see Equations ( f5Ti) . ( [55]) ). 

We deal with the evolution of constraints ( 168]) and (169]) . in the following section. 

6. f(R) Constraint propagation equations in pure 3+1 language 

Note that Equations ( |76]) and ( |77j) involve both spacetime and purely spatial objects. 
This is not a form that easily yields a comparison between the constraint propagation 
equations obtained in the Frittelli approach with those obtained in the 3+1 approach. 
Nor is it convenient for integration in a 3+1 numerical implementation that could serve 
as a check of the numerical integration of the evolution equations of the dynamical 
variables. For this reason, we now cast these equations in pure 3+1 language. To our 
knowledge such a calculation has never been published before, hence it is instructive to 
include it here. 

§ We note here that Eq. (11) in 85 , differs from our Equation ([73]) by a factor of 2 in the term 
Hn^V ^n a . We believe that this discrepancy is simply due to a typographical error. 
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Alternative expressions for the extrinsic curvature are (see e.g. Equations (2.49), 
(2.52) in [78]) 

D( a np) = -K a/3 and K aP = -V a n^ - n a ap, (78) 

where a a = n^V pn a is the acceleration of normal observers, also equal to (see Eq. (2.22) 
in [93]) 

a/3 = Dp In a. (79) 
From Equation (1781) it can be shown that 

D p n p = -K and D p n a = -K p a . (80) 
By use of Equations ( 1781) . ( 1791) and ( IHUl) . Equations ( 1761) and ( ITTj) can be written as 
n^V^H = - D^M" + (1 + X)HK 

-2M u D u lna, (81) 
n^V^M u = - im^D^H + n u M^D^ In a + ikTK 

+ M^/-(l + A)ii\D' / lna. (82) 

The identities, V Q if = <9 Q iJ, ^ U D^H = j^d^H, M a D a lna = M a d a \na, V M ikP = 

can be used to replace the covariant derivatives that occur 
above. Furthermore, the timelike unit vector (n M ) can be replaced by Equation ( 16"T1) . 
Equations (IHTj) and (182|) can then be written as 

d t H = ftdiH - 2M i d i a - aDM + (1 + fi)aHK, (83) 

d t M j = - ^ ij diH + ftdiMi - (4) r^ M i 

+ {4) Ti k M i /3 k + n j M l dia + aM j K 

+ aM l Ki - (1 + fi)^ ij Hdia, (84) 

where we have focused on the spatial indices of M M , since M M is purely spatial. 

Using DiM 1 = Y^diMj + 'ijM^, and the expressions for the Lie derivatives of 
the constraints along an^ 

£ an H =d t H-f3%H, (85) 
£ an M j = d t M j - ftdiM 3 + Wdifi (86) 

we write Equations (155!) and ([541) as 



£ Qn iJ = _ 2M i d i a - af j diMj 

+ a^V k i5 M k + (1 + n)otHK, (87) 

£ an M j = - X^diH + A j i M i + aikFif 
-(l + A) 7 u '^a + 7 jfc M^/3 fc 

- M E ^ m d asm , (88) 
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where 

A j i = - (4) r| + ^Y\ k P k - a-^dta + aK(. (89) 

We now need to express the Christoffel symbols associated with the spacetime 
metric g pu , that appear in Equation (1591 . in terms of the 3-metric 7^ and the gauge 
variables. We do this as follows 

(4) r^ = 2 gJP ( dig °P + do9i P ~ 9 P gi0 "> 

= ^9 3 °digoo + ^'(digat + d g ie - d £ g i0 ). (90) 



Using the relations between g^ v and a, (3\ 7^ 

g Q0 = -a? + p l p e , <?ck = A, ^ = 7^-a"W, (91) 
Equation (|90|) eventually becomes 

- i( a - 2 /W3i7* - y'^ft - 7 J '^o7tf) 

- \{l j %Pi + a- 2 f3 j f3%7u - a~W^A), (92) 
where we have also used the following identities 

^ = 7% d k y j = -i is i jm d klsm . (93) 

The next object that appears in Equation (I8"9"j) . and which we cast in 3+1 language 
is P k ^T^ ik . We can write this as 

[3 k ^P lk = \p k g ip {d i9kp + d k g ip - d p g ik ) 

= ^ k g J °(dig k0 + d k g i0 - d g lk ) 

+ \(3 k g ] \d igM + d k g u - d e g ik ) . (94) 
By virtue of Equations (!9T|) and f|93|) . Equation fl94"|) finally becomes 

-±a- 2 p k p i d " fik + p k r* ik 

- a- 2 pip k p e r' ik , (95) 



or equivalently 



P kii) T\ k = l -a' 2 p k ^d i p k + ±a~ 2 p k pid k p i 
- l - a -ip k pidv lik + p k V ik 



1 



W^Tm, (96) 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 



18 



where T^ ik stand for the Christoffel symbols associated with the 3-metric. 
By use of Equations fl92|) and fl96|) . Equation (189]) becomes 

A j i = - - + 

+ f3 k P ik + olK\. (97) 
From the evolution equation of the 3-metric, Equation (|60|) . we have 

= -aK\ + -y'^ft + ~7^A - 7^^. (98) 

Substitution of Equation f )98|) into Equation (19?]) yields 

A i. = + ^ e (3 k d aik + 2aKh- (99) 

Finally, substituting Equation ( 19"9~|) into Equation ([8"8"]) yields the desired result, 

£ an M j = - A7 ij '^i7 + 2aK{M l + aM j K 

- (1 + X)^ ij Hd ia . (100) 

Equations flH7|) and (11 001) are the hamiltonian and momentum constraint 
propagation equations in pure 3+1 language, where the Lie derivatives are given in 
Equations ( I85JI and flggp . 

We have already established that the form of the constraint propagation equations 
is the same for both vacuum and non-vacuum spacetimes, and is independent of the 
form of the function f(R). Thus, to validate our equations we can use known results 
that apply to the case of GR, and have been derived by using the "brute force" method. 

For this reason we now compare our results with results published in [831 EH] that 
apply for f(R) = R, i.e., for the Einstein equations. In these two papers the evolution 
equations of the constraints were presented assuming = 0. In [83] the 3+1 approach 
was employed to derive the constraint propagation equations. For direct comparison 
with these published results we also derive the evolution equations for Ti = 2H and the 
evolution for Mj which were used in [831 EE] instead. 

Using 

£ an Mi = M J £ an jij + ^ij£ an M 3 

= - 2aK ij M 1 + lij £ an M j (101) 

and replacing H = %/2 in ( 1H71) and f llOOp we obtain the following alternative form for 
the constraint propagation equations 

£ an n = - Wdia - 2a^ ij d i M j 

+ 2a^T k ij M k + (1 + X)aHK, (102) 

£ an Mi = —XdiU + aMiK - (1 + \)-Hdia. (103) 

For A = 1 Equations (jlU2j) and (I103p become precisely the same as the expressions 
in [83], when the quantities, C ki j = dkjij — D ki j, defined in [83], satisfy C k ij = 0. In 
that work C k ij are constraints that arise from the introduction of the auxiliary variables 
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Df-ij = dk'jij, which were used to reduce the ADM formulation to lst-order. Also, 
a straightforward calculation shows that the expressions above are equivalent to the 
corresponding expressions in [86]. From Equations (I102j) and (I103p it is again evident 
that the constraints remain satisfied (Ji = Mi = 0), if they are initially satisfied. 

We now turn our attention to the evolution equations of Ci and C#. To derive 
the evolution of Cj and Cy we simply take a time derivative of Ci and C^, use the 
commutation relation did t = d t di, and replace the time derivative of variables via 
Equations ( 164]) and ( 166]) to find that 

dtCt = (3 s C st (104) 

and 

dtCij = £/sCij, 

where 

l- ,(•,, = niCj + ( ■,!<); r + Cji-r. (106) 

Equations (11041) and (11051) imply that if the constraints Ci, C^ are initially satisfied, 
then the evolution equations will preserve the constraints. 

We stress again that Equations (11021) . (I103p . (I104p and (11051) are valid not only for 
vacuum, but also for non- vacuum spacetimes, as well as for any viable f(R) model. This 
is a new result that to our knowledge has not been pointed out previously and is not 
trivial to prove, if one employs the 3+1 or "brute force" method to derive the constraint 
propagation equations. Here we proved this without prior knowledge of the evolution 
equations of the dynamical variables Kij,jij, on the basis of the Frittelli approach. 

Finally, we note that the agreement between our expressions and results obtained 
by the "brute force" approach confirms that Equation ( 1731) is correct (see footnote in 
page [15]) . 



(105) 



7. Constraint Propagation in the Jordan and Einstein frames 

The form of the f{R) field equations both in the Jordan frame (see Equations ffTUl) 
and ( TT2|) ) and in the Einstein frame (see Equations ( |2T]) and ( 123]) ) is the same as the 
BD formulation of f{R) gravity. For this reason, it is evident from our discussion 
in Section [5] that the form of the constraint propagation equations in these two 
frames must be the same as those in the BD formulation, provided that we define 
the hamiltonian and momentum constraints analogously to Equations fl5Tl) and (155]) 
with one important caveat; The foliation in the Einstein frame must be based on the 
conformal (Einstein) metric and not the physical (Jordan), i.e., the induced 3- metric 
on spacelike hypersurfaces must be j a b = g a b + n a fib, where the normal timelike vector 
now satisfies gab'n a n b — — 1 and not g a bn a n h = —1. Note that this last condition is 
only a mathematical requirement for the 3+1 machinery to remain the same. Physical 
conclusions must still be drawn based on the Jordan metric. 
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Finally, we note that if one applies the general recipe for a 3+1 decomposition (see 
Section H]) to more general scalar-tensor theories of gravity considered in [82], then the 
constraint propagation equations will be the same as our Equations (11021) . (11031) . (11041) 
and (11051) . This is because the covariant (Jordan frame) formulation of these theories 
obtains the same form as Equations ( 1321) and (1341) that we considered here, which lead 
to same decomposition fl59l) . 

8. Summary and Discussion 

We have extended the ADM constraint propagation equations, using the Frittelli method 
[85J, to generic metric f(R) gravity represented as a BD theory. For direct comparison 
with published results, we wrote our general evolution equations of the constraints 
(defined via Equations ( !5lT) and ( |55l) ) in the same form as the original equations given 
in [85]. This mathematical form, given by Equations (JTBl) and ( 1771) . combines both 
spacetime and purely spatial objects. To make transparent the connection between these 
equations and the language of the 3+1 decomposition of spacetime, we cast Equations 
(176]) and (1771) in pure 3+1 form, i.e., in a form that involves only scalar and purely 
spatial objects and their derivatives (see Equations (11021) and (1103!) ). The 3+1 form is 
the mathematical form the evolution equations of the constraints would take on, if one 
had employed a "brute force" 3+1 approach for performing this derivation. The brute 
force approach requires prior knowledge of the exact 3+1 equations and is much more 
involved. 

The main result of this work is that the mathematical form of the constraint 
propagation equations is the same for both vacuum and non-vacuum spacetimes, as 
well as for any (viable) form of the function f{R), provided that Tffl and T$y (see 
Section Hj) are absorbed properly in the definition of the constraints. We have also 
argued that the mathematical form of the evolution of the constraints for 3+1 f(R) 
gravity in the Jordan frame remains the same as that of the BD-equivalent 3+1 f{R) 
gravity. This result holds true in the Einstein frame, too, if the spacetime foliation 
is chosen based on the Einstein metric g^ v and not the physical (Jordan) metric g^ v . 
Finally, a comparison between our equations and previous GR results, using the 3+1 
approach, shows that all expressions for the constraint propagation equations agree. 

We end this work by pointing out that the 3+1 BD-equivalent f(R) equations can 
be incorporated in current numerical relativity codes with only minor additional effort. 
For example, the minimum requirement for studying vacuum spacetimes in viable f(R) 
models is to include the contribution of the scalar field stress-energy tensor T$ (see 
Section H]) and implement a scalar field solver in the form of Equations ( !64l) -( l66l) . As 
in GR, it is almost certain that the stability of numerical implementations of the fully 
non-linear equations of f{R) gravity will be sensitive to the formulation used. Given 
that the structure of the f(R) constraint propagation equations is fundamentally the 
same as that of the ADM formulation, we believe that dynamical f(R) simulations will 
benefit from formulations such as the Baumgarte-Shapiro-Shibata-Nakamura (BSSN; 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 21 



[HUES]) approach or the generalized harmonic decomposition [HE]- If these formalisms 
fail, other approaches such as those proposed in [93, [HU [98] may prove useful. We hope 
that this work will serve as a starting point for relativists to develop fully dynamical 
codes for viable f(R) models. 

Acknowledgments 

We would like to thank Carlos Cunha for useful conversations. This paper was supported 
in part by NSF Grants PHY06-50377 and PHY09-63136 as well as NASA Grants 
NNX07AG96G and NNX10AI73G to the University of Illinois at Urbana-Champaign. 
Ignacy Sawicki acknowledges support by the DFG through TRR33 "The Dark Universe" . 

References 

[1] Clifford M. Will. The confrontation between general relativity and experiment. Living Reviews 
in Relativity, 9(3), 2006. 

[2] A. G. Riess et al. Observational evidence from supernovae for an accelerating universe and a 

cosmological constant. Astron. J., 116:1009-1038, 1998. 
[3] S. Perlmutter et al. Measurements of ft and A from 42 High-Redshift Supernovae. Astrophys. J., 

517:565-586, 1999. 

[4] P. J. E. Peebles and B. Ratra. The cosmological constant and dark energy. Rev. Mod. Phys., 
75:559-606, 2003. 

[5] Sean M. Carroll. The cosmological constant. Living Reviews in Relativity, 4(1), 2001. 

[6] E. Komatsu et al. Seven- Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: 

Cosmological Interpretation. ArXiv e-prints, astro-ph/1001.4538, 2010. 
[7] R. Amanullah et al. Spectra and Light Curves of Six Type la Supernovae at 0.511 < z < 1.12 

and the Union2 Compilation. Ap. J., 716:712-738, 2010. 
[8] Beth A. Reid et al. Cosmological Constraints from the Clustering of the Sloan Digital Sky Survey 

DR7 Luminous Red Galaxies. Mon. Not. Roy. Astron. Soc., 404:60-85, 2010. 
[9] Will J. Percival et al. Baryon Acoustic Oscillations in the Sloan Digital Sky Survey Data Release 
7 Galaxy Sample. Mon. Not. Roy. Astron. Soc., 401:2148-2168, 2010. 
[10] Edmund J. Copeland, M. Sami, and Shinji Tsujikawa. Dynamics of dark energy. Int. J. Mod. 

Phys., D15T753-1936, 2006. 
[11] Antonio De Felice and Shinji Tsujikawa. f(R) Theories. Living Reviews in Relativity, 13(3), 2010. 
[12] T. P. Sotiriou and V. Faraoni. f(R) theories of gravity. Rev. Mod. Phys., 82:451-497, 2010. 
[13] Philip D. Mannheim. Alternatives to dark matter and dark energy. Progress in Particle and 

Nuclear Physics, 56(2) :340 - 445, 2006. 
[14] Roy Maartens and Kazuya Koyama. Brane-world gravity. Living Reviews in Relativity, 13(5), 
2010. 

[15] Luca Amendola, Kari Enqvist, and Tomi Koivisto. Unifying Einstein and Palatini gravities. ArXiv 

e-prints, gr-qc/1010.4776, 2010. 
[16] P. G. Bergmann. Comments on the scalar-tensor theory. Int. J. Theor. Phys., 1:25-36, 1968. 
[17] T. V. Ruzmaikina and A. A. Ruzmaikin. Quadratic corrections to the lagrangian density of the 

gravitational field and the singularity. Zh. Eksp. Teor. Fiz., 57:680, 1969. 
[18] H. A. Buchdahl. Non-linear lagrangians and cosmological theory. Mon. Not. R. Astron. Soc, 

150:1-8, 1970. 

[19] A. A. Starobinsky. A new type of isotropic cosmological models without singularity. Phys. Lett. 
B, 91:99-102, 1980. 

[20] S. Capozziello. Curvature quintessence. Int. J. Mod. Phys. D, 11:483-491, 2002. 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 22 

[21] S. Capozzicllo, S. Carloni, and A. Troisi. Quintessence without scalar fields. In Recent Research 
Developments in Astronomy and Astrophysics 1, page 625. Research Signpost, Trivandrum, 
India, 2003. 

[22] S. Capozzicllo, V. F. Cardone, S. Carloni, and A. Troisi. Curvature quintessence matched with 

observational data. Int. J. Mod. Phys. D, 12:1969-1982, 2003. 
[23] S. M. Carroll, V. Duvvuri, M. Trodden, and M. S. Turner. Is cosmic speed-up due to new 

gravitational physics? Phys. Rev. D, 70:043528, 2004. 
[24] S. Nojiri and S. D. Odintsov. Modified gravity with negative and positive powers of the curvature: 

Unification of the inflation and of the cosmic acceleration. Phys. Rev. D, 68:123512, 2003. 
[25] T. Chiba. 1/R gravity and scalar-tensor gravity. Phys. Lett. B, 575:1-3, 2003. 
[26] G. J. Olmo. The gravity lagrangian according to solar system experiments. Phys. Rev. Lett., 

95:261102, 2005. 

[27] G. J. Olmo. Post-Newtonian constraints on f(R) cosmologies in metric and Palatini formalism. 

Phys. Rev. D, 72:083505, 2005. 
[28] L. Amcndola, D. Polarski, and S. Tsujikawa. Are f(R) dark energy models cosmologically viable? 

Phys. Rev. Lett., 98:131302, 2007. 
[29] L. Amendola, D. Polarski, and S. Tsujikawa. Power-laws f(R) theories are cosmologically 

unacceptable. Int. J. Mod. Phys. D, 16:1555-1561, 2007. 
[30] L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa. Conditions for the cosmological viability 

of f(R) dark energy models. Phys. Rev. D, 75:083504, 2007. 
[31] Thomas Faulkner, Max Tegmark, Emory F. Bunn, and Yi Mao. Constraining f(R) gravity as a 

scalar tensor theory. Phys. Rev., D76:063505, 2007. 
[32] W. Hu and I. Sawicki. Models of f(R) Cosmic Acceleration that Evade Solar-System Tests. Phys. 

Rev. D, 76:064004, 2007. 

[33] Hiroaki Oyaizu. Non-linear evolution of f(R) cosmologies I: methodology. Phys. Rev., D78T23523, 
2008. 

[34] Hiroaki Oyaizu, Marcos Vinicius Lima, and Wayne Hu. Non-linear evolution of f(R) cosmologies 

II: power spectrum. Phys. Rev., D78T23524, 2008. 
[35] Fabian Schmidt, Marcos Vinicius Lima, Hiroaki Oyaizu, and Wayne Hu. Non-linear Evolution of 

f(R) Cosmologies III: Halo Statistics. Phys. Rev., D79:083518, 2009. 
[36] Simone Ferraro, Fabian Schmidt, and Wayne Hu. Cluster Abundance in f(R) Gravity Models. 

ArXiv e-pmnts, astro-ph/1011.0992, 2010. 
[37] Gong-Bo Zhao, Baojiu Li, and Kazuya Koyama. N-body Simulations for f(R) Gravity using a 

Self-adaptive Particle-Mesh Code. ArXiv e-prints, astro-ph/1011.1257, 2010. 
[38] B. Li and J. D. Barrow. The Cosmology of f(R) Gravity in the Metric Variational Approach. 

Phys. Rev. D, 75:084010, 2007. 
[39] A. A. Starobinsky. Disappearing cosmological constant in f(R) gravity. J. Exp. Theor. Phys. 

Lett, 86:157-163, 2007. 

[40] S. A. Appleby and R. A. Battye. Do consistent f(R) models mimic general relativity plus A? 

Phys. Lett. B, 654:7-12, 2007. 
[41] N. Deruelle, M. Sasaki, and Y. Sendouda. 'Detuned' f(R) gravity and dark energy. Phys. Rev. 

D, 77:124024, 2008. 

[42] G. Cognola et al. A class of viable modified f(R) gravities describing inflation and the onset of 

accelerated expansion. Phys. Rev. D, 77:046009, 2008. 
[43] E. V. Lindcr. Exponential gravity. Phys. Rev. D, 80:123528, 2009. 

[44] A. V. Frolov. A Singularity Problem with f(R) Dark Energy. Phys. Rev. Lett, 101:061103, 2008. 
[45] T. Kobayashi and K. Maeda. Relativistic stars in f(R) gravity, and absence thereof. Phys. Rev. 
D, 78:064019, 2008. 

[46] E. Babichev and D. Langlois. Relativistic stars in f(R) gravity. Phys. Rev. D, 80:121501, 2009. 
[47] A. Upadhye and W. Hu. The existence of relativistic stars in f(R) gravity. Phys. Rev. D, 
80:064002, 2009. 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 23 



[48] E. Babichev and D. Langlois. Relativistic stars in f(R) and scalar-tensor theories. Phys. Rev. D, 
81:124051, 2010. 

[49] D. Lai, F. A. Rasio, and S. L. Shapiro. Ellipsoidal figures of equilibrium - Compressible models. 
ApJS, 88:205-252, 1993. 

[50] D. Lai, F. A. Rasio, and S. L. Shapiro. Hydrodynamic instability and coalescence of binary neutron 

stars. Ap. J., 420:811-829, 1994. 
[51] S. Chandrasekhar. Ellipsoidal figures of equilibrium. The Silliman Foundation Lectures, New 

Haven: Yale University Press, 1969, 1969. 
[52] M. Shibata, T. W. Baumgarte, and S. L. Shapiro. The Bar-Mode Instability in Differentially 

Rotating Neutron Stars: Simulations in Full General Relativity. Ap. J., 542:453-463, 2000. 
[53] M. Saijo, M. Shibata, T. W. Baumgarte, and S. L. Shapiro. Dynamical Bar Instability in Rotating 

Stars: Effect of General Relativity. Ap. J., 548:919-931, 2001. 
[54] Kostas D. Kokkotas and Bernd Schmidt. Quasi- normal modes of stars and black holes. Living 

Reviews in Relativity, 2(2), 1999. 
[55] M. D. Duez. Numerical relativity confronts compact neutron star binaries: a review and status 

report. Classical and Quantum Gravity, 27(11):114002 — h, June 2010. 
[56] I. Hinder. The current status of binary black hole simulations in numerical relativity. Classical 

and Quantum Gravity, 27(11):114004-+, June 2010. 
[57] E. Rantsiou, S. Kobayashi, P. Laguna, and F. A. Rasio. Mergers of Black Hole-Neutron Star 

Binaries. I. Methods and First Results. Ap. J., 680:1326-1349, 2008. 
[58] F. Loffler, L. Rezzolla, and M. Ansorg. Numerical evolutions of a black hole-neutron star system 

in full general relativity: Head-on collision. Phys. Rev. D, 74(10): 104018 ■+, 2006. 
[59] J. A. Faber, T. W. Baumgarte, S. L. Shapiro, K. Taniguchi, and F. A. Rasio. Dynamical evolution 

of black hole-neutron star binaries in general relativity: Simulations of tidal disruption. Phys. 

Rev. D, 73(2):024012 +, 2006. 
[60] J. A. Faber, T. W. Baumgarte, S. L. Shapiro, and K. Taniguchi. General Relativistic Binary 

Merger Simulations and Short Gamma-Ray Bursts. Ap. J. Letters, 64LL93-L96, 2006. 
[61] M. Shibata and K. Uryii. Merger of black hole-neutron star binaries: Nonspinning black hole case. 

Phys. Rev. D, 74(12):121503-+, 2006. 
[62] M. Shibata and K. Uryu. Merger of black hole neutron star binaries in full general relativity. 

Classical and Quantum Gravity, 24:125 — h, 2007. 
[63] Masaru Shibata and Keisuke Taniguchi. Merger of black hole and neutron star in general relativity: 

Tidal disruption, torus mass, and gravitational waves. Phys. Rev. D, 77(8):084015, 2008. 
[64] T. Yamamoto, M. Shibata, and K. Taniguchi. Simulating coalescing compact binaries by a new 

code (SACRA). Phys. Rev. D, 78(6):064054-+, 2008. 
[65] Z. B. Etiennc, J. A. Faber, Y. T. Liu, S. L. Shapiro, K. Taniguchi, and T. W. Baumgarte. Fully 

general relativistic simulations of black hole-neutron star mergers. Phys. Rev. D, 77(8):084002- 

+, 2008. 

[66] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baumgarte. General relativistic simulations 
of black-hole-neutron-star mergers: Effects of black-hole spin. Phys. Rev. D, 79(4):044024 — h, 
2009. 

[67] M. D. Duez et al. Evolving black hole-neutron star binaries in general relativity using 
pseudospectral and finite difference methods. Phys. Rev. D, 78(10):104015 — h, November 2008. 

[68] M. Shibata, K. Kyutoku, T. Yamamoto, and K. Taniguchi. Gravitational waves from black hole- 
neutron star binaries: Classification of waveforms. Phys. Rev. D, 79(4):044030 — h, February 
2009. 

[69] K. Kyutoku, M. Shibata, and K. Taniguchi. Quasiequilibrium states of black hole-neutron star 
binaries in the moving-puncture framework. Phys. Rev. D, 79(12):124018 — h, June 2009. 

[70] P. M. Motl et al. Fully Relativistic Simulations of the Inspiral and Merger of Black Hole - Neutron 
Star Binaries. In Bulletin of the American Astronomical Society, volume 41 of Bulletin of the 
American Astronomical Society, pages 295 — h, January 2010. 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 24 



[71] S. Chawla et al. Mergers of Magnetized Neutron Stars with Spinning Black Holes: Disruption, 

Accretion and Fallback. ArXiv e-prints, gr-qc/1006.2839, 2010. 
[72] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, and S. A. Teukolsky. Equation of state effects in 

black hole-neutron star mergers. Classical and Quantum Gravity, 27(11):114106 — h, June 2010. 
[73] F. Pannarale, A. Tonita, and L. Rczzolla. Black hole-neutron star mergers and short GRBs: a 

relativistic toy model to estimate the mass of the torus. ArXiv e-prints, astro-ph/1007.4160, 

2010. 

[74] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukolsky. Black hole-neutron star mergers: 
effects of the orientation of the black hole spin. ArXiv e-prints, astro-ph/1007.4203, 2010. 

[75] K. Kyutoku, M. Shibata, and K. Taniguchi. Gravitational waves from nonspinning black hole- 
neutron star binaries: dependence on equations of state. ArXiv e-prints, astro-ph/1008.1460, 
2010. 

[76] Vasileios Paschalidis, Morgan MacLeod, Thomas W. Baumgarte, and Stuart L. Shapiro. Merger 

of white dwarf-neutron star binaries: Prelude to hydrodynamic simulations in general relativity. 

Phys. Rev. D, 80(2):024006, 2009. 
[77] V. Paschalidis, Z. Etienne, Y. T. Liu, and S. L. Shapiro. Head-on collisions of binary white 

dwarf-neutron stars: Simulations in full general relativity. ArXiv e-prints, astro-ph/1009.4932, 

2010. 

[78] Thomas W. Baumgarte and S. L. Shapiro. Numerical Relativity: Solving Einstein's Equations on 

the Computer. Cambridge University Press, 2010. 
[79] Miguel Alcubierre. Introduction to 3+1 Numerical Relativity. Oxford University Press, 2008. 
[80] C. Brans and R. H. Dicke. Mach's principle and a relativistic theory of gravitation. Phys. Rev., 

124:925-935, 1961. 

[81] Marcelo Salgado. The cauchy problem of scalar-tensor theories of gravity. Classical and Quantum 

Gravity, 23(14) :4719, 2006. 
[82] N. Lanahan-Tremblay and V. Faraoni. The Cauchy problem of f(R) gravity. Classical and 

Quantum Gravity, 24:5667-5679, November 2007. 
[83] Vasileios Paschalidis. Mixed hyperbolic - second-order-parabolic formulations of general relativity. 

Phys. Rev. D, 78:024002, 2008. 
[84] Gioel Calabrese. A remedy for constraint growth in numerical relativity: the Maxwell case. 

Classical and Quantum Gravity, 21(17):4025, 2004. 
[85] S. Frittelli. Note on the propagation of the constraints in standard 3+1 general relativity. Phys. 

Rev. D, 55:5992-5996, 1997. 
[86] Gen Yoneda and Hisa-aki Shinkai. Constraint propagation in the family of ADM systems. Class. 

and Quant. Grav., 19:1027, 2002. 
[87] H.-a. Shinkai and G. Yoneda. Letter: Constraint Propagation in (N + l)-Dimensional Space- Time. 

General Relativity and Gravitation, 36:1931-1937, August 2004. 
[88] Steven Weinberg. Gravitation and Cosmology: Principles and Applications of the General Theory 

of Relativity. John Wiley & Sons, Inc., 1972. 
[89] R. H. Dicke. Mach's principle and invariance under transformation of units. Phys. Rev., 125:2163- 

2167, 1962. 

[90] N. Deruelle, Y. Sendouda, and A. Youssef. Various Hamiltonian formulations of f(R) gravity and 
their canonical relationships. Phys. Rev. D, 80(8):084032-+, October 2009. 

[91] R. Arnowitt, S. Deser, and C.W. Misner. Gravitation: An Introduction to Current Research. 
Cambridge University Press, Wiley, New York, 1962. 

[92] J. W. York. Sources of Gravitational Radiation. Cambridge University Press, Cambridge, England, 
1979. 

[93] T. W. Baumgarte and S. L. Shapiro. Numerical relativity and compact binaries. Phys. Rep., 
376:41-131, 2003. 

[94] Masaru Shibata and Takashi Nakamura. Evolution of three-dimensional gravitational waves: 
Harmonic slicing case. Phys. Rev. D, 52(10) :5428-5444, Nov 1995. 



Constraint propagation equations of the 3+1 decomposition of f(R) gravity 25 

[95] Thomas W. Baumgarte and Stuart L. Shapiro. Numerical integration of Einstein's field equations. 

Phys. Rev. D, 59(2):024007, Dec 1998. 
[96] F. Pretorius. Evolution of Binary Black- Hole Spacetimes. Physical Review Letters, 95(12):121101- 

+, 2005. 

[97] Vasileios Paschalidis, Jakob Hansen, and Alexei Khokhlov. Numerical performance of the 
parabolized ADM formulation of general relativity. Phys. Rev. D, 78(6):064048, Sep 2008. 

[98] David R. Fiske. Toward making the constraint hypersurface an attractor in free evolution. Phys. 
Rev. D, 69(4):047501, Feb 2004. 



